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Abstract. Stress patterns in static granular media exhibit unusual features 
when compared to either liquids or elastic solids. Qualitatively, we attribute 
these features to the presence of 'stress paths', whose geometry depends on 
the construction history and controls the propagation of stresses. Stress 
paths can cause random focussing of stresses (large fluctuations) as well as 
systematic deflections (arching). We describe simple physical models that 
capture some of these effects. In these models, the 'stress paths' become 
identified with the characteristic 'light rays' of wavelike (hyperbolic) equa- 
tions for force propagation. Such models account for the 'pressure dip' below 
conical sandpiles built by pouring from a point source, and explain qual- 
itatively the large stress fluctuations observed experimentally in granular 
matter. The differences between this approach and more conventional mod- 
elling strategies (based on elastoplastic or rigid-plastic models) are high- 
lighted, focusing on the role of boundary conditions. Our models provide a 
continuum picture in which granular materials are viewed as fragile mat- 
ter, able to support without rearranging only a subset of the static external 
loadings admissible for a normal elastic solid. 
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1. Introduction 

Stress patterns in granular media exhibit some rather unusual features when 
compared to either liquids or elastic solids. For example, the vertical pres- 
sure below conical sandpiles does not follow the height of material above 
a particular point, but rather has a minimum underneath the apex of the 
pile [1, 2, 3]. Furthermore, local stress fluctuations are large, sometimes 
on length scales much larger than the grain size. For example, repeatedly 
pouring the very same amount of powder in a silo results in fluctuations 
of the weight supported by the bottom plate of 20% or more [4, 5]. Quali- 
tatively, these features are attributed to the presence of stress paths which 
can focus the stress field into localized regions and also deflect it to cause 
"arching" (see also [6, 7], and [8] for early qualitative experiments). 

More quantitative experiments were recently performed by Liu et al. [9], 
Brockbank et al. [2] and Mueth et al. [10], where the local fluctuations of 
the normal stress deep inside a silo or at the base of a sandpile were mea- 
sured. It was found that the stress probability distribution is rather broad, 
decaying exponentially for large stresses. A simple 'scalar' (one component) 
model for stress propagation was introduced and studied in detail [9, 11], 
which predicts a stress probability distribution in good agreement with ex- 
perimental (and numerical) data. However, this model only considers the 
vertical normal component of the stress tensor, and discards all the other 
components. 

A fully 'tensorial' model for stress propagation in homogeneous granular 
media was proposed in [12, 3, 13] to account for the pressure 'dip' described 
above. The most striking feature of this model is that the stress propagation 
equation is (at least in two dimensions) a wave equation, with the vertical 
axis playing the role of time. In this model, 'stress paths' naturally appear 
as the characteristics - or the 'light rays' - of the corresponding hyperbolic 
equation. Note that the standard equations of elasticity are elliptic; the 
fundamental difference between these two cases will be discussed later. Both 
must also be contrasted with the scalar model, which corresponds to a 
parabolic equation, and in which stresses travel almost vertically. 

The aim of this paper is to review some of the recent theoretical work in 
this field. We will start by summarizing the content of the scalar '<?-modeP, 
which, although unsatisfactory in several respects, offers the advantage of 
simplicity. Keeping the spirit of the 'g-model', we then show that the in- 
troduction of the shear stress fundamentally modifies the structure of the 
equations and leads to wave-like propagation. Several issues concerning the 
solution of this wave equation are then discussed. We then consider some 
variants of the wave equation, in particular to account for local inhomo- 
geneities or for anisotropy, for example induced by the construction history. 
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(The importance of construction history in granular matter has been ac- 
knowledged for at least a century [14].) Among a family of models, the 
'Fixed Principal Axis' (fpa) limit plays a special role which we discuss in 
relation with experimental data on conical sandpiles. We then consider in 
more detail the differences between this approach and some more conven- 
tional modelling strategies (based on elastoplastic or rigid-plastic models), 
and try to shed some light on the controversy that the "hyperbolic" ap- 
proach to sandpile modelling has raised recently [15, 16, 17]. This discussion 
focuses in particular on the role of boundary conditions. 

Much of our discussion will, for clarity, be limited to two dimensional 
piles. However, most of the recent experimental results are in three di- 
mensions; indeed it was in this context that the physical assumptions of 
Refs.[3, 13] were made. 



2. The Scalar Model 

The main assumption of the scalar model is that only the vertical nor- 
mal component of the stress tensor w = o zz (the 'weight') needs to be 
considered. Supposing for simplicity that the grains reside on the nodes 
of a two-dimensional lattice (see figure 1), the simplest model for weight 
propagation is: 

= w g + q+(i-l,j-l)w{i-l,j-l) + q-(i + l,j-l)w(i+l,j-l) (1) 

where l w g ^ is the weight of a grain, and q±(i, j) are 'transmission' coefficients 
giving the fraction of weight which the grain transmits to its right 

(resp. left) neighbour immediately below. Mass conservation imposes that 
Q+(hj) + Q-(hj) = 1 f° r an h j' s - The case of an ordered pile of identical 
grains would correspond to q± = \. In this case, the equation (1) describes 
the 'time'-evolution (along the j direction) of the probability density for 
a random walker on a line. The authors of [9, 11] then propose to take 
into account the local disorder in packing, grain sizes and shapes, etc., 
by choosing q + (i,j) to be independent random numbers (subject to the 
above constraint), for example uniformly distributed between and 1. This 
case is interesting because it leads to an exact solution for the local weight 
distribution P(w), in the limit j — > oo: 

i» = ^e*p-- (2) 

where 2W = jw g is the average weight. Liu et al. [9, 11] have argued that 
the exponential tail for large w is generic; however, if the maximum value 
of q is qu < 1, one can show that P(w) decays faster than an exponential: 

log PH OC^oo -J* = (3) 

log 2q M 
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Figure 1. The 'q-model' model with two neighbours. q±'s are independent random 
variables, satisfying the weight conservation constraint: q+(i,j) + q-(i,j) = 1. 



(Notice that (3 = 1 whenever qu = 1, and that (3 = oo in the ordered case 
qu = 1/2)- In this sense the exponential tail of P*(w) is not universal: it 
requires the possibility that one of the q can be arbitrarily close to 1, i.e. 
that there is a non zero probability that one grain is entirely bearing on 
one of its downward neighbours (local arching). 

In the continuum limit, the q- model is equivalent to a diffusion equation 
with a random convection term (related to q + — (?_), for which several results 
are known. However, on large scales, the random convection term merely 
renormalises the diffusion constant; hence, in this model, the response to 
a localised overload at the top of a pile spreads out diffusively as V DH, 
where H is the height of the pile and D is the diffusion constant, which 
is of the order of the grain size a. In the limit H » a, the spreading is 
negligible, which means that the weight essentially propagates vertically. 
Hence, the g-model predicts a 'hump' in the pressure profile underneath a 
sandpile, directly reflecting its shape. 

A way to accomodate a "pressure dip" within this scalar picture was 
suggested by Edwards (although in a slightly different language). Consider 
for example a sandpile built from a point source: the history of the grains 
will certainly inprint a certain oriented 'texture' to the contact network, 
which can be modelled, within the present scalar model, as a nonzero mean 
value of the 'convection term', the sign of which depends on which side 
of the pile is chosen. Let us call Vq the average value of this term on the 
x > side of the pile, with — Vq on the other side. The differential equation 
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describing propagation now reads, in the absence of disorder: 

d t w + d x [V sign(x)w] = p + D d xx w (4) 

Solving this equation in a sandpile geometry leads to a weight minimum 
around x = 0. Equation (4) gives a precise mathematical content to Ed- 
wards' idea of arching in sandpiles [18], as the physical mechanism leading 
to a 'dip' in the pressure distribution [1]. As discussed below (see also 
[3, 13]), this can be taken much further within a tensorial framework. Let 
us also note that (4) can in fact be obtained naturally within an extended 
g-model model, with an extra rule accounting for the fact that a grain 
can slide and lose contact with one of its two downward neighbours [19]. 
However, this extra sliding rule implicitly refers to the existence of shear 
stresses, absent in the scalar model. It is more satisfactory to recognize from 
the outset that stress has a tensorial, rather than scalar, nature. Models 
which do this are decribed next. 

3. A Tensorial Model 

3.1. THE WAVE EQUATION 

It is useful to start with a simple toy model for stress propagation, which 
is the analogue of the model presented in figure 1. We now consider the 
case of three downward neighbours (see figure 2), for a reason which will 
become clear below. Each grain transmits to its downward neighbours not 




Figure 2. Three neighbour configuration. Each grain tranmits two force components to 
its downward neighbours. A fraction p of the vertical component is transmitted through 
the middle leg. 

one, but two force components: one along the vertical axis t and one along 
x, which we call respectively F t (i,j) and F x (i,j). For simplicity, we assume 
that each 'leg' emerging from a given grain can only transport the force 
parallel to itself (but more general rules could be invented). Assuming that 
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the transmission rules are locally symmetric, and that a fraction p < 1 of 
the vertical component travels through the middle leg, we find: 

F x (i,j) = ^[F x (i-l,j-l)+F x (i + l,j-l)] 

+ 1(1 _ p ) tanV [^(t - 1, J - 1) - F t (i + 1)} (5) 

F t (i,j) = w +pF t (i,j-l) + ^(l-p)[F t (i-l,j-l)+F t (i + l,j-l)] 

+ + (6) 

where ip is the angle between grains, defined in figure 2. Taking the contin- 
uum limit of the above equations leads to: 

8 t F t + d x F x = p (7) 
d t F x + d x [clF t ] = (8) 

where Cq = (1 — p) tan 2 ip . Eliminating (say) F x between the above two 
equations leads to a wave equation for Ft, where the vertical coordinate t 
plays the role of time and cq is the equivalent of the 'speed of light'. In 
particular, the stress does not propagate vertically, as it does in the scalar 
model, but rather along two rays, each at a non zero angle ±(p such that 
Co = tan cp. Note that <p ^ ip in general (unless p = 0); the angle at which 
stress propagates has nothing to do with the underlying lattice structure 
and can take any value depending on the local rules for force transmission. 
We chose a three- leg model to illustrate this particular point. 

The above derivation can be reformulated in terms of classical contin- 
uum mechanics as follows. Considering all stress tensor components o~ij, the 
equilibrium equation reads, 

d t o- t t + d x o xt = p (9) 
d t o- tx + d x u xx = (10) 

Identifying the local average of Ft with o~tt and that of F x with o~t x , we see 
that the above equations (7, 8) are actually identical to (9, 10) provided 
o~tx = o~ x t (which corresponds to the absence of local torque) and 

0- XX = Cq o t t (11) 

This relation between normal stresses was postulated in [12] as the sim- 
plest "constitutive relation" among stress components, obeying the correct 
symmetries, that one can possibly assume. The term "constitutive relation" 
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normally refers to a relation between stress and strain, but the model under 
discussion has no strain variables defined; instead the particles are viewed 
as completely rigid. (Equations (9, 10) are then indeterminate unless a fur- 
ther hypothesis relating the stresses themselves is made.) This particular 
choice can be interpreted as a local Janssen approximation [20] . We return 
later to a more detailed discussion of closure equations of this type. In the 
present case, the parameter Cq must encode relevant details of the local 
geometry of the packing (friction, shape of grains, etc.) and may thereby 
depend on the construction history of the grain assembly. Only for simple, 
'homogeneous' histories (such as constructing a uniform sandbed using a 
sieve) will c\ be everywhere constant on the mesoscopic scale. Even then, 
unless an ordered packing is somehow created, local fluctuations of Cg will 
always be present. 

3.2. SOME SIMPLE SITUATIONS 

The simplest situation is that of an infinitely wide layer of sand, of depth 
H, with a localized (<5-function) overload at the top. The weight at the 
bottom then defines the response function of the wave equation, which, in 
two dimensions, is the sum of two 5 peaks localised at x = ±coH. 

Next, one can consider the sandpile geometry. For a pile at repose, the 
position of the free surfaces are x = ±cz, where c = cot <fi with (f> the repose 
angle. On these surfaces, all the stresses vanish. This boundary condition 
is then (for given cq and c) sufficient to solve for the stress field everywhere 
in the pile. (See Section 8 below.) One then find that the vertical normal 
component of the stress is piecewise linear as a function of x. In particular, 
for —cqH < x < coH, o t t is constant. Therefore, in two dimensions, this 
model [12] predicts a flat-topped stress profile rather than a dip. 

For a pile created by depositing grains from above (for example by 
sieving sand onto a disc) it is natural to expect the free surface to be a 
slip plane. (This is a plane across which the stress components saturate the 
Mohr-Coulomb condition.) Interestingly, this provides a relation between 
Co and the friction angle 4>, which reads: Cg = 1/(1 + 2tan 2 <^>) (note that 
since c = 1/ tan0, one has automatically c > Co). Under these conditions 
one finds that the 'plastic' region (where the Mohr-Coulomb condition is 
saturated) extend inward from the surface to encompass the outer 'wings' 
of the pile (i.e. cqz < |x| < cz); see figure 3. This follows from the solution 
of the model and is not an a priori assumption, of the kind commonly 
made in elastoplastic modelling (e.g., [16]). In three dimensions, a second 
closure relation is required [12], but in all cases the stress profile has a 
broad maximum at the center of the pile. Now, however, the Mohr-Coulomb 
condition is only saturated in the immediate vicinity of the free surface - 
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the 'plastic' region has zero volume in three dimensions [12, 13]. 
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Figure 3. Reduced normal (upper) and shear stress (lower) curves against S = r tan <j>/t 
for the symmetric wave equation ("BCC") and the fpa model in two dimensions. Inset: 
the yield function Y. The Mohr- Coulomb inequality is saturated when Y = 1 



4. Symmetries and Constitutive Relations 

Although above it was motivated in the context of a specific microscopic 
model, the linear constitutive relation (11) can be viewed, independently of 
any microscopic model, as the simplest closure equation compatible with the 
symmetries of the problem. The latter include a local reflection symmetry in 
which x — xq is changed to xq — x (with xq an arbitrary reflection plane) and 
also a form of "dilational" symmetry known as RSF ("radial stress field") 
scaling. RSF scaling depends on the absence of any characteristic stress 
scale, which follows if the Young's modulus of the grains is sufficiently much 
larger than any stresses arising in the granular assembly being studied. Such 
scaling, which requires the stress distributions beneath piles of different 
heights to have the same shape, is quite well confirmed in some (but not 
all) experiments on conical sandpiles [1, 2, 15]. 

Even with these two symmetries, one can consider more complicated 
(nonlinear) constitutive relations among stresses, which must be of the 
form [12]: 



Note that the Mohr-Coulomb condition itself can be written in this form. 
Viewed as a constitutive equation, it defines a rigid-plastic model whose 




(12) 
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physical content is to assume that, everywhere in the material, a plane can 
be found across which slip failure is about to occur. (The name "incipient 
failure everywhere", ife, aptly describes this model [12, 3, 13].) All closures 
of the form (12) lead to a hyperbolic equations for stresses, although in the 
general case the characteristic directions of propagation (the 'light rays' of 
the corresponding wave equation) depend on the loading and therefore vary 
with position. 

An interesting situation arises when local reflection symmetry is broken. 
This is the case, for example, in sandpiles created by pouring from a point 
source onto a rough surface - which is the usual mode of construction. In 
such a pile, all grains arriving at the apex of the pile roll (in two dimen- 
sions) either to the right or to the left. The two halves of the pile therefore 
have different construction histories that are mirror images of each other. 
This violates local reflection symmetry, and in general permits constitutive 
equations such as: 

a xx = cla a g(^^) (13) 

The simplest case (found e.g. by expanding Q to first order in the shear to 
normal stress ratio) corresponds to a family of (quasi-) linear constitutive 
relations [13]: 

<?xx = cl<j tt + H sign(x)cr xt (14) 

The previous, symmetrical, case has n = 0. For nonzero fi, (14) again leads 
to a wave equation, although this time anisotropic, in the sense that the 
two characteristic rays make asymmetric angles to the vertical axis. Note 
that such a model can be obtained from rules such as those in figure 2 sim- 
ply by having an asymmetric partitioning of forces between the supporting 
grains (or indeed by tilting the entire packing). Note also that x = is a 
singular line across which the directions of propagation change discontinu- 
ously. Microscopically, /x ^ also leads to an unequal sharing of the weight 
of a grain between the two characteristic rays propagating downward from 
it. For jjL < 0, most of the weight travel outwards; this provides, within a 
fully tensorial model, a mathematical description of the tendency to form 
arches, as developed by Edwards for the scalar case. 

Solving these anisotropic wave equations for sandpiles in two dimensions 
one again finds for a tt a piecewise linear function, which now has a sharp 
maximum at x = when /j, > 0, but a minimum for fj, < 0, in accord with 
the arching scenario mentioned above (see figure 3). If one furthermore 
imposes, as above, that the free surfaces are slip planes, one finds a relation 
between Cq, \i and <p. 



C ° = l + 2tan^ [1 - Man ^ ] 



(15) 
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One again finds the result that the material throughout the outer wings 
of the pile (exterior to the triangle formed by the characteristics passing 
through the apex) are at incipient (Mohr-Coulomb) failure. 

4.1. THE FPA MODEL 

Among possible forms of (14) there is a particular case, corresponding to 
Cq = 1, which has some intriguing special properties: 

&xx = ott ~ 2sign(x) tan(0) a xt (16) 

(The resulting value of ^ is that appropriate to the boundary conditions 
stated above.) Specifically in this case, the principal axis of the stress tensor 
coincide with the 'light rays' and hence have a fixed direction throughout 
the pile (up to a reflection symmetry across the line x = 0) . This particular 
case was called the fpa ("fixed principal axes") constitutive equation in 
[3, 13]. It corresponds to an assumption that the anisotropic texture of the 
medium imparts a local stress rule which requires that the orientation of 
the stress tensor is fixed at the time of burial and, thereafter, cannot change 
under further loading (The values of the stress components themselves can, 
of course, change.). Since all material elements are buried near the free 
surface of the pile, where the stress tensor orientation is known the fact that 
the surface is itself a slip plane [3] , this fixes the principal axes throughout 
the pile. The resulting stresses are shown in figure 3. 

At first sight, there is a problem with the FPA description on the sym- 
metry axis of the pile (x = 0), where the stress ellipsoid is required si- 
multaneously to have two conflicting orientations. However, since Cq = 1 
the stress is isotropic at the centre, and there is no conflict. However, this 
shows that the defining feature of the FPA model is very easily lost; as soon 
as a slightly different Cg is chosen, the principal axes are no longer of fixed 
orientation but rotate smoothly as one passes from one side of the pile to 
the other. 

5. Experiments on Cones and Wedges 

We do not have space here for an extended discussion of the experimental 
data on cones and wedges; this can be found elsewhere [21]. Here we merely 
summarize some of the most important points. 

5.1. CONES 

The extension of the fpa model to three dimensions, like the other mod- 
els discussed above, requires a second constitutive relation among stresses 
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to close the problem. Several well-known candidates [22] exist for this sec- 
ondary closure and give similar results; these results are in surprisingly 
good agreement with the experimental data for the pressure dip in three 
dimensional conical sandpiles built by pouring from a point source onto a 
rough rigid support [3, 13, 2]; see figure 4. 
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Figure 4- Comparison of fpa model (using a uniaxial secondary closure [12, 3, 13]) with 
scaled experimental data of Smid and Novosad and (*) that of Brockbank et al (averaged 
over three piles). Upper and lower curves denote normal and shear stresses. The data 
is used to calculate the total weight of the pile which is then used as a scale factor for 
stresses. The horizontal coordinate is scaled by the pile radius. 

We believe that the experimental data strongly suggest an FPA-like 
model (corresponding to Cq ~ 1, Eq.16), although they do not prove that 
the principal axes of a material element are necessarily fixed. Of greater 
physical significance is the idea that the characteristic rays of the wave 
equation could be fixed at burial. (Below we will connect this with the idea 
of stress paths.) This is the content of a more general interpretation (the 
"oriented stress linearity" model of [13]) in which \i (or Cq) is a parameter 
characterising the local anisotropy of the granular packing. This is fixed by 
construction history, in a manner as yet unspecified; but for a point-source 
pile values close to the FPA limit are suggested by the data. In contrast, for 
a pile constructed by sieving, we might expect the local reflection symme- 
try to be unbroken {p, = 0), leading instead to a smooth maximum in the 
pressure, as predicted from Eq.ll. Experiments on the construction-history 
dependence of the stress profile would therefore be welcome. 

Note that the idea that some property of the medium, represented by 
a constitutive equation among stresses, is "fixed at burial" is an important 
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simplifying assumption of our approach to sandpile modelling. This has 
been called the "perfect memory" assumption [13]. 



5.2. WEDGES: SURFACE AVALANCHES OR DEEP YIELD? 

Savage [15] has pointed out that some classical data on wedges, as opposed 
to cones, are apparently at odds with these ideas. Such wedges are of tri- 
angular cross-section but very long in the third direction, i.e., quasi-two 
dimensional. According to classical reports [23, 24], the stress distribution 
beneath such wedges shows little dependence on the construction history, 
but a very strong dependence on whether the base supporting the wedge is 
allowed to sag. Without basal sag, almost no pressure dip is reported [15]. 

At first sight, then, the pressure dip in conical piles might also be at- 
tributable to basal sag. We strongly believe that this viewpoint is unten- 
able, especially in the light of the data of Brockbank et al [2] (see figure 
4) in which the measurement system leads to discernable localized indenta- 
tion, rather than a slight curvature (sagging) of the base under regions of 
high pressure. Moreover, the classical data on wedges is somewhat scant, 
of dubious accuracy, and in most cases does not specify the construction 
history of the wedge in any detail. We would therefore encourage renewed 
experimental investigation of wedges of sand. 

There is, in any case, an important difference in the physics of wedges 
and cones. This concerns the geometry of the "plastic" region (that in which 
the Mohr-Coulomb inequality is saturated) near the surface of the pile. All 
our models, including the FPA limit [13], predict that this is infinitely narrow 
in the three dimensional conical pile, but that for a wedge it extends through 
a large outer zone (exterior to the two characteristic rays passing through 
the apex). In the first case, the perfect memory assumption appears self- 
consistent: the presence of a thin yield layer at the surface suggests that the 
pile grows by surface avalanches which do not disturb its internal structure 
too much. In contrast, the status of the perfect memory assumption is, 
for a wedge, far less clear. Since a broad zone of marginal instability exists 
beneath the surface, the surface flow could cause "deep yield" events [25, 17] 
which would disrupt the internal structure of the pile [21]. This could lead 
to a local isotropisation of the granular texture and cause values of fi in (14) 
closer to those of the symmetric propagation model (Eq.ll) than the fpa 
model (Eq.16). Much more detailed experimental data is certainly needed, 
however, before these ideas can be put to the test. 
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6. The Role of Local Inhomogeneities 

6.1. A STOCHASTIC WAVE EQUATION 

Provided that local conservation laws (those arising from mechanical equi- 
librium) are obeyed, many local rules for force transmission are a priori 
compatible with the existence of contacts among rigid particles [26, 27]. 
Therefore, even if there is a definite mean relationship among stresses at 
the meso-scale (as models like fpa assume), one can expect randomness in 
the local transmission coefficients. The simplest model for this and other 
sources of randomness is to introduce a randomly varying 'speed of light' 
Co- This could describe the fact that, for example, the parameter p in the 
model of figure 2 can vary from grain to grain. 

This suggests the following stochastic wave equation for stress propa- 
gation in two dimensions: 



where v(x, t) is a random noise. We have studied this equation in great 
details [28] using perturbation theory, and find that for weak disorder, the 
average response function now has two peaks of finite (diffusive) width 
(rather than two 5 peaks in the zero disorder case), and that the 'speed 
of light' is renormalized to a lower value cr. More interestingly, the un- 
averaged response function takes negative (and rather large) values. This 
may be of crucial importance since it suggests a fundamental instability 
of granular matter to external perturbations. Suppose indeed that as a re- 
sult of a distant perturbation, a certain grain receives a negative (upward) 
force larger than the pre-existing downward vertical pressure. This grain 
will then move and a local rearrangement of contacts will occur. If stability 
is to be recovered, this rearrangement must induce a variation of co(x, t) so 
as to reduce the cause of the instability. Thus, the stochastic wave equa- 
tion implicitly demands rules similar to those introduced in [19] to describe 
extreme sensitivity to external perturbations in silos. The present model, 
which is purely static, does not say what happens when a local rearrange- 
ment occurs, but certainly suggests that small perturbations will induce 
such rearrangements. One can show that typically a perturbation of order 
the weight of one grain is enough to oblige rearrangements somewhere else 
in the pile [28]. This would cease to be true if large enough overload was 
applied to ensure that that all grains are subject to a vertical compression 
greatly in excess of their weight. But it is not clear that such an overload 
is ever really possible: even at the bottom of a deep pile, a finite fraction 
of grains may be effectively non-loadbearing [10]. 

We have also studied the weight-weight correlation function, in the ge- 
ometry of a flat layer with a random overload at the top, neglecting gravity. 



dtto-tt = dxx Cq(1 + v(x, t)) a tt 



(17) 
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We find [28] that, as a function of (horizontal) distance, the correlation has 
two peaks. The first one is of course at separation Ax = 0, while the second 
is at Ax = 2crH, which simply means that two points at the bottom of 
the packing connected by 'light rays' to the same point on the top, share 
information about the overload. This result is of importance since the shape 
of this correlation function clearly differs from the corresponding one in the 
scalar model, which is a single peak at Ax = 0. It also differs from that 
pertaining to a simple elastic medium, where the weight-weight correlation, 
in this geometry, decays only very slowly: correlations extend to the scale of 
the system size itself. Measuring carefully the averaged correlation function 
of a granular system under an overload could then confirm (or disprove) 
the presence of a ray-like propagation. A stress correlation function was 
recently measured in [10] and found to be featureless, but measurements 
extended only to very short lengthscales: Ax < 5a, as compared to the 
height of the pile H ~ 100a. We thus expect the features of the correlation 
function to show up on much larger scales (~ 2crH) than those measured 
so far. 



6.2. STRESS HISTOGRAM 



We have seen above that within a scalar approach, an exponential-like 
distribution (possibly of the type exp — w^, with (3 > 1) is expected [9, 
11]. One can ask whether this exponential distribution survives within a 
tensorial description. So far, we only have partial numerical results, based 
on a direct simulation of the three-leg model introduced above, with a 
random p chosen between and pm- This scheme is thus very close in 
spirit to the g-model. However, as emphasized above, the local vertical 
forces are not everywhere positive; one should thus introduce an extra rule 
to cope with this instability. Several possibilities come to mind, but we have 
not yet explored them (see however [26, 27]). Nevertheless, the large force 
region, which is presumably not sensitive to the presence of negative forces, 
behaves much in the same way as in the scalar g-model. In particular, the 
tail of the distribution decays as exp— w@, with [3 ~ 1 when pu = 1, and 
with (3 > 1 when pm < 1- More work is needed to understand the physical 
implications of the presence of negative forces and any relation this may 
have to the static avalanche phenomenon [19]. However, the above results 
show that the tail of the force distribution is only exponential in a 'strong 
disorder' limit, where local 'arching' (i.e. one grain entirely bearing on a 
single downward neighbour) has a nonzero probability. 
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7. Modelling Strategies for Static Granular Media 

We now compare our approach (as outlined above) to the problem of sand- 
pile modelling with previous ones, and address various criticisms of it that 
have recently been made. Specifically we are interested in calculating the 
stress distribution at the base of a pile constructed on a rough, rigid sup- 
port. For clarity of the discussion, we again limit the mathematics to two 
dimensions, although we emphasize again that our work, particularly that 
on the FPA model, was developed in the context of three dimensional conical 
piles. 

7.1. LOCAL RULES FOR FORCE TRANSFER 

Above we have described an approach to the modelling of stress propaga- 
tion in granular media based on local rules for the transfer of forces between 
grains. Broadly speaking, to leading order in a gradient expansion [21, 28] 
the closure (14) exhausts the possibilities for local models in which the 
forces passed from a grain to its neighbours in the layer below involve a 
linear decomposition of the "incident force" (F x ,Ft), which is taken as the 
vector sum of forces acting on the grain from those in the layer above. 
Somewhat similar considerations underlie the so-called "clastic" theories of 
"discontinua" introduced by Trollope [29]. Indeed, some authors have con- 
fused the two models, and our approach has been criticized for "reinventing 
theories... that are well-known in another disciplinary area" [15]. 

However, Trollope's model, though linear, does not lead to Eq.14. Its 
distinguishing feature is instead that the vector sum of the incident forces 
on a grain is not taken before applying a rule to determine the outgo- 
ing forces from that grain. The outgoing forces instead depend separately 
on each of the incident force contributions. This feature is, in our view, 
strongly unphysical: if a grain is subjected to two small extra forces, whose 
vector sum is vertical, from its neighbours in the layer above, the force 
increments exerted on the grains below should be equivalent to a small in- 
crease in its weight. Within Trollope's rules, this is not the case [30, 21]. It is 
therefore disingenuous to criticize our approach on the grounds that 'Trol- 
lope's clastic or discontinua model, that was rejected by Wittmer et al. on 
the grounds of being "unphysical", actually contains the fpa solution' [15]. 
For, although Trollope's model can be tuned to give the same stress pattern 
as the fpa model in a symmetric two dimensional pile, these two models 
have distinctly different physical content and give differing predictions in 
other geometries. 
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7.2. HYPERBOLIC CONTINUUM MODELS 

Although motivated in part by simple packing models (e.g. figure 2), our 
approach does not require such a specific microscopic interpretation and 
can instead by formulated in purely continuum mechanical terms [31]. This 
is done by postulating a constitutive relations among stresses such as (14). 
(The latter includes, for special values of fj,, both (11) and (16) which de- 
scribe respectively the case of symmetric stress propagation and the fpa 
model.) We have noted already the existence of an alternative constitutive 
equation among stresses, corresponding to the Mohr-Coulomb rigid-plastic 
model (or ife model) which is a widely studied continuum theory of gran- 
ular media. In fact this reads (in two dimensions): 



All these approaches lead to hyperbolic equations for stress propagation. 

A reasonable question [15] is to ask why we are not satisfied with the 
ife closure (18). The main cause is that we can see no physical reason for 
it to be true. Indeed, even supporters of this rigid-plastic model do not 
usually propose it as an accurate description of sandpile behaviour; it is 
more often viewed as a way of generating certain "limit-state" solutions. 
In the simplest geometries these solutions correspond to taking the — or 
+ sign in (18); these are often referred to as "active" and "passive" limit- 
states, respectively [32]. It is easily established [13, 21] that for a sandpile 
at its repose angle, only one solution of the resulting equations exists in 
which the sign choice is everywhere the same. This is the active solution, 
and it shows a hump, not a dip, in the vertical normal stress beneath the 
apex. Savage, however, draws attention to a "passive" solution, having a 
pronounced dip beneath the apex [15]. This solution actually contains a 
pair of matching planes between an inner region where the positive root of 
(18) is taken, and an outer region where the negative is chosen. 

In principle this does not exhaust the repertoire of ife solutions for the 
sandpile: there should exist others, with larger numbers of matching planes 
between segments of alternating sign [21]. In any case, the predictive power 
of the rigid-plastic approach largely depends on a belief that the limit-state 
solutions can be "generally regarded as bounds between which other states 
can exist, i.e., when the material is behaving in an elastic or elastoplas- 
tic manner" [15]. Unfortunately, this belief is unfounded: counterexamples 
exist (figure 5), even among elastoplastic models of the simplest kind [21]. 
Therefore the so-called limit-states should be regarded merely as "rule of 
thumb" estimates; these may be useful for engineering purposes, but do not 
shed much light on the physics of stress propagation in granular matter. 




(18) 
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Figure 5. Vertical normal stress found from Eq.ll (the BCC model, [12]), for a pile at 
angle of repose <f> — 30 degrees, compared to the active and passive ife solutions. (The 
ife solutions are obtained by shooting from the midplane for P = ia tt + a xx )/2 and 
the polar angle 9 as functions of the direction of the principal axis *&.) Note that active 
and passive ife solutions do not bound the stress, either in the model of Eq.ll or in the 
clastoplastic model of [16], which, for a certain parameter choice, yields identical results 
to Eq.ll (see [16]). 



7.3. ELASTOPLASTIC MODELS 

All the models considered above make no mention of strain variables. A 
partial justification for this was given in Refs.[3, 13], namely that the ex- 
perimental data obey radial stress-field (RSF) scaling, which implies that 
there is no characteristic length-scale. Since elastic deformation under grav- 
ity introduces such a length-scale (a "sagging length") the observation of 
RSF scaling to experimental accuracy in most but not all the data [1, 2, 15] 
suggests that elastic deformation is not significant. 

This does not imply that an elastic or elastoplastic description of sand- 
piles is impossible; but it shows that in any such description, the limit of 
a large elastic modulus appears to be the relevant one. This limit yields 
equations, in the bulk of the medium, for which strain variables cancel out; 
and this fact is usually exploited in elastoplastic calculations (see, e.g., [16]). 
Note that it is tempting, but entirely wrong, to assume that strain variables 
on the boundary of the medium also cancel in this limit (see below). 

In fact, as correctly noted by Savage [15] results similar to those of the 
fpa model and the related models described above can, in two dimension 
(only), be obtained within such an elastoplastic modelling approach. Typ- 
ically, an inner, linear elastic region is matched, by hand, onto an outer 
plastic one. An example of this procedure was described recently by Can- 
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telaube and Goddard [16] whose approach is similar to earlier work by 
Samsioe [33]. This analysis can be made to give mathematically identical 
results to those found with some, if not all, values of \x in Eq.14. 

As they stand, however, such elastoplastic analyses are devoid of phys- 
ical meaning, for the following reason. Recall that the aim of the exercise 
is to calculate the forces measured at the base of a pile of sand. Recall 
also the well-known theorem that to find the equilibrium state of an elastic 
body, one must specify either the surface force field or the displacement 
field at all points on its boundary [34]. Accordingly, it is meaningless to 
"calculate" the forces at the base of an elastoplastic pile without specifying 
a boundary condition at the bottom surface of any elastic zones present. 
To specify as boundary conditions the forces themselves there is mathe- 
matically legitimate, but cannot really be described as a "solution" to the 
sandpile problem. On the other hand, to specify the displacements is (by 
this stage) formally impossible, since the relevant variables have already 
been eliminated by taking the limit of a high elastic modulus. 

The only way to circumvent this difficulty is to specify the displacement 
field at the base first, and then take the limit of a high modulus afterwards. 
To obtain finite forces, the displacements must be allowed to tend to zero 
but, crucially, the results depend on how this limit is taken. This is illus- 
trated in figure 6. The challenge to elastoplastic modellers is then to decide 
what (infinitesimal) displacement field should be chosen at the base. For a 
genuine elastoplastic body, which is placed on a rough rigid support, this 
displacement field depends on the precise manner in which the body was 
brought into its state of rest. For a sandpile constructed by (say) pouring 
sand grainwise from a point source the problem seems much less well-posed 
since (despite assertions to the contrary [15]) there is no obvious physical 
definition of the displacement field in such a pile [35]. 

7.4. ELASTIC INDETERMINACY 

The circumstances just described apply to any model of a sandpile in which 
the elliptic equations of elasticity are invoked in all or part of the medium. 
We have called this the problem of "elastic indeterminacy" [21]. Several 
responses to this challenge are possible. One is to assume that the elasto- 
plastic description is essentially sound, and seek some physical procedure 
whereby the displacement field at the support, or some equivalent informa- 
tion, is determined [36]. This is certainly worth pursuing, but we suspect 
that any description which does not explicitly consider the construction 
history of the pile is very unlikely to be successful. 

A second response, which appears to be that of Evesque [17] is to con- 
clude that the experimental results are and must be indeterminate. Put 
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Figure 6. Consider the static equilibrium of an elastic body of finite modulus resting 
on a completely rough surface. Starting from any initial configuration, another can be 
generated by pulling and pushing parts of the body horizontally across the base {i.e., 
changing the displacements there); if this is rough, the new state will still be in static 
equilibrium. This will generate a stress distribution, across the supporting surface and 
within the pile, that differs from the original one. If the limit of a large modulus is now 
taken (at fixed stress), this procedure allows one to generate arbitrary differences in the 
stress distribution while generating neither finite distortions in the shape of the body, 
nor any forces at its free surface. An exactly analagous "elastic indeterminacy" exists in 
any simple elastoplastic theory of sandpiles, where an elastic zone, in contact with part 
of the base, is attached at matching surfaces to a plastic zone. 



differently, Evesque argues that the external forces acting on the base of 
a pile can in fact be varied at will by the experimentalist. Certainly, if 
sandpiles obey Hooke's law, he must be right: for an elastic body, the ex- 
perimentalist is free to prod about at the base of a pile in such a way as to 
change arbitrarily the forces acting there. (Moreover, in the large modulus 
limit, only infinitesimal proddings are required.) 

In contrast to this, experimental reports clearly suggest [1, 2] that the 
forces on the base can be measured more or less reproducibly, and (though 
subject to statistical fluctuations) do not vary too much from one pile to 
another. Moreover, the experimental data for the forces show RSF scaling at 
the base; this is confirmed by comparing data from piles of different heights. 
Within an elastoplastic framework, this scaling should be violated by the 
(arbitrary) boundary forces at the base and hence can only be expected in 
the upper extremity of the pile [16, 21]. 

All this suggests a third response to the challenge of elastic indeter- 
minacy. This is to argue that Hooke's law has very little relevance to the 
mechanics of sandpiles and that models based on local force propagation 
rules among grains are far closer to the real physics of the problem. The 
models we have developed along these lines give hyperbolic equations for 
stress propagation, and in doing so contradict Hooke's law in a fundamental 
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way. This applies even in their incremental response to small added loads. 
Whether this is a drawback or an advantage depends on one's view of the 
physics. Certainly, if one believes that any granular assembly must behave 
elastically under sufficiently small incremental loads, then models such as 
ours can only describe the behaviour beyond some finite threshold [36]. 
(We are as yet unconvinced of whether this threshold is finite; for exam- 
ple it might vanish in the high modulus limit.) In any case, in calculating 
the response to gravity itself we believe that any such threshold is easily 
exceeded throughout the pile, and hence that our approach has far more 
to offer than models invoking elastic deformations from some hypothetical 
unstrained (i.e., zero gravity) state. On the other hand, the existence of a 
threshold might make it somewhat harder to detect experimentally the re- 
sponse and correlation functions [12, 13, 28] which are, as described above, 
strong signatures of hyperbolic stress propagation laws. 

The varying perspectives laid out above are not necessarily completely 
contradictory, in that the global ('coarse-grained') features of the stress 
pattern could be governed by determinate equations but the details not. 
As well as being subject to elastic indeterminacy, the latter could be affected 
by randomness in the grain packings, which may be exquisitely sensitive to 
temperature and other poorly-controlled parameters [19]. 

8. Boundary Conditions in Hyperbolic Models 

We now consider more carefully the role of boundary conditions in our 
approach. In contrast to the elliptic equations of elasticity, the hyperbolic 
equations arising from constitutive relations among stresses admit definite 
solutions for the stresses at the base of a freestanding pile. By the same 
token, if one tries to apply all the boundary conditions appropriate to an 
elastic body (e.g. specifying the surface forces acting over the entire sur- 
face), then in general no solution will exist of the hyperbolic equations that 
pertain to any particular choice of constitutive relation. We discuss these 
two aspects in turn. 

8.1. DETERMINACY OF THE SANDPILE 

The procedure for a sandpile is of course to specify zero-force boundary con- 
ditions at the free (upper) surfaces. Within our description, the response 
arising from a localized body force (a "source term" in the equations) prop- 
agates downward along two characteristics passing through the source. In 
models obeying (14) these characteristics are, in addition, straight lines. 
The force on the base is found simply by summing the contributions from 
all the body forces; this is a fully determinate procedure for any closed 
set of hyperbolic equations [12, 3, 13]. Note that in principle, one could 
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envisage propagation also along the "backward" characteristics (see figure 
7(a)). This is forbidden since these cut the free surface; any such propaga- 
tion can only arise in the presence of a nonzero surface force, in violation of 
the boundary conditions. Therefore the fact that propagation occurs only 
along downward characteristics is not related to the fact that gravity acts 
downward; it arises because we know what forces act at the free surface (the 
forces there are zero). Suppose we had instead an inverse problem: a pile 
or bed with some unspecified overload at the upper surface, for which the 
forces acting at the base had been measured. In this case, the information 
from the known forces could be propagated along the upward characteristics 
to find the unknown overload. 

(a) (b) V 




A 

Figure 7. (a) The response to a localized force is found by resolving it along characteris- 
tics through the point of application, propagating along those which do not cut a surface 
on which the relevant force component is specified. For a pile under gravity, propagation 
is only along the downward rays, (b) Admissible boundary conditions cannot specify sep- 
arately the force component at both ends of the same characteristic. If these forces are 
unbalanced (after allowing for body forces), static equilibrium is lost. 



8.2. STRESS PATHS 

Informally speaking, the hyperbolic problem is determined once half of the 
boundary forces are specified. More precisely (figure 7(b)) one is required 
to specify the surface force tangential to each characteristic ray, at one 
end and one end only. The corresponding force acting at the other end 
is obliged to balance it, allowing for any body forces acting tangentially 
along the ray. If it does not do so, then within our modelling approach, 
the material ceases to be in static equilibrium. This is no different from the 
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corresponding statement for a fluid or liquid crystal; if boundary conditions 
are applied that violate the conditions for static equilibrium, some sort of 
motion results. Unlike a fluid, however, for a granular medium we expect 
such motion to be in the form of a finite rearrangement rather than a steady 
flow. Such a rearrangement will change the microtexture of the material, 
and thereby alter the constitutive relation among stresses. We expect it to 
do so in such a way that the new constitutive relation is compatible with 
the imposed forces. 

Although simplified, we believe that this picture correctly captures some 
of the essential physics of stress paths. Such paths are load-bearing struc- 
tures within the contact network and, in the simplest approximation of 
straight, unbranched paths, these must have the property described above 
(the forces on two ends of a path must balance) . Stress paths should there- 
fore be identified (on the average) with the characteristics of our hyperbolic 
equations. The mean orientation of the stress paths is then reflected in a 
constitutive equation such as (11) or (16). 

The physics of our modelling approach is thus to assume that the mean 
orientation of stress paths, in each element the material, is fixed at burial. 
(This does not necessarily require that the individual paths are themselves 
fixed.) This was called earlier the "perfect memory assumption". We think 
it reasonable to assume that the stress paths will not change their average 
orientations so long as they are able to support the applied load. But if a 
load is applied which they cannot support, rearrangement is inevitable [17]. 
This causes some part of the pile, at least, to adopt a new microtexture and 
thereby a new constitutive relation. In other words, loadings of this kind 
cause the construction history of the pile to change. 

9. Conclusion: Sandpiles as Fragile Matter 

Within our modelling approach, a granular assembly is able to support 
some, but not all, of the surface loads that would be supportable by an elas- 
tic continuum. Such models may therefore provide an interesting paradigm 
for the behaviour of "fragile matter", a concept which may be useful in 
other systems where certain combinations of applied forces, even if small, 
are enough to force irreversible reconstruction of the material. (Such sys- 
tems could include a number of disordered soft solid materials such as defect 
textures in liquid crystals.) 

In the present context, fragility arises from the the requirement of tan- 
gential force balance along stress paths. If this is violated at the boundary, 
even infinitesimally, then internal rearrangement must occur, causing new 
stress paths to form, so as to support the load [37]. Obviously, this might 
be rather too simple a picture - for example, branching of stress paths is 
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ignored. Thus it remains possible that Hookean behaviour is recovered for 
sufficiently small perturbing forces [36, 17]. However, for practical purposes 
we believe our approach, by capturing at least some of the physics of stress 
paths, may have rather more to offer than ideas based on the physics of 
conventional (homogenous, isotropic) elastic, elastoplastic or rigid-plastic 
media. Phenomena that seem to be addressable in such terms include that 
of arching and, in its various aspects, the role of noise. 

Clearly, there is much scope for developing these models further. In par- 
ticular it would be very useful to have an understanding of the crossover, 
if one indeed exists, between fragile and elastic regimes: the latter should 
ultimately be restored in a sufficiently large pile (beyond the "sagging" 
length). Equally important is to confront these and other modelling ap- 
proaches with much more demanding experimental tests, of which several 
were suggested above. 
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